clear;

listDir = dir('Data*');
saveArg = 1;
numBins = 72;
intrudeErupt = zeros(length(listDir),1);
time = 183:50:583;
MadMedTheshold = 0.6;


%% Processing

L = nan(length(listDir),1);
B = L;
thetaMean = L;
thetaMedian = L;
thetaDev = L;
thetaMad = L;
correctedMean = L;
correctedMedian = L;
correctedDev = L;
correctedMad = L;
reducedMean = L;
reducedMedian = L;
reducedDev = L;
reducedMad = L;
reducedCorMean = L;
reducedCorMedian = L;
reducedCorDev = L;
reducedCorMad = L;
vectorDistMed = nan(numBins,length(listDir));
vectorDistMad = nan(numBins,length(listDir));
correctedDistMed = nan(numBins,length(listDir));
correctedDistMad = nan(numBins,length(listDir));
reducedDistMed = nan(numBins,length(listDir));
reducedDistMad = nan(numBins,length(listDir));
reducedCorrDistMed = nan(numBins,length(listDir));
reducedCorrDistMad = nan(numBins,length(listDir));
velocityMean = L;

for i = 1:length(listDir)
    data = importdata(sprintf('%s\\B00001.dat',listDir(i).name));
    data = data.data;
    x = data(:,1);
    y = data(:,2);
    vx = data(:,3);
    vy = data(:,4);
    
    noInd = find(vx==0&vy==0);
    x(noInd) = [];
    y(noInd) = [];
    vx(noInd) = [];
    vy(noInd) = [];
    L(i) = max(y) - min(y);
    B(i) = max(x) - min(x);
    
    vector = sqrt(vx.^2+vy.^2);
    theta = acos(vx./vector);
    
    
    threshold = max([10^-5,median(vector)-mad(vector)]);
    
    xG = linspace(min(x),max(x),70);
    yG = linspace(min(y),max(y),70);
    UG = nan(length(yG),length(xG));
    VG = UG;
    devUG = UG;
    devVG = VG;
    for i1 = 1:length(yG)-1
    for i2 = 1:length(xG)-1
        checkInd = x>=xG(i2) & x<=xG(i2+1) & ... 
            y>=yG(i1) & y<=yG(i1+1) & vector>threshold;
        UG(i1,i2) = median(vx(checkInd));
        VG(i1,i2) = median(vy(checkInd));
        devUG(i1,i2) = mad(vx(checkInd));
        devVG(i1,i2) = mad(vy(checkInd));
    end
    end
%     figure(1000); clf;
%     subplot(1,2,1);
%     pcolor(log10(abs(devUG)./abs(UG))); caxis([0,1]);
%     subplot(1,2,2);
%     pcolor(log10(abs(devVG)./abs(VG))); caxis([0,1]);
    
    UG(abs(devUG)./abs(UG)>MadMedTheshold) = NaN;
    VG(abs(devVG)./abs(VG)>MadMedTheshold) = NaN;
    
%     figure(1000); clf;
%     subplot(1,2,1);
%     pcolor(log10(abs(devUG)./abs(UG))); caxis([0,1]);
%     subplot(1,2,2);
%     pcolor(log10(abs(devVG)./abs(VG))); caxis([0,1]);
    
    
    reducedVector = sqrt(UG.^2+VG.^2);
    reducedTheta = acos(UG./reducedVector);
    
    correctedTheta = theta;
    correctedTheta(vy<0) = 2*pi - correctedTheta(vy<0);
    correctedTheta = correctedTheta + pi/2;
    correctedTheta(correctedTheta>2*pi) = ...
        correctedTheta(correctedTheta>2*pi) - 2*pi;
    
    reducedCorrected = reducedTheta;
    reducedCorrected(VG<0) = 2*pi - reducedCorrected(VG<0);
    reducedCorrected = reducedCorrected + pi/2;
    reducedCorrected(reducedCorrected>2*pi) = ...
        reducedCorrected(reducedCorrected>2*pi) - 2*pi;
    
    
    for j = 2:numBins
        checkInd = find(theta>=2*pi*(j-1)/numBins & ...
           theta<=2*pi*j/numBins);
        vectorDistMed(j,i) = median(vector(checkInd));
        vectorDistMad(j,i) = mad(vector(checkInd));
        
        checkInd = find(correctedTheta>=2*pi*(j-1)/numBins & ...
           correctedTheta<=2*pi*j/numBins);
        correctedDistMed(j,i) = median(vector(checkInd));
        correctedDistMad(j,i) = mad(vector(checkInd));
        
        checkInd = find(reducedTheta>=2*pi*(j-1)/numBins & ...
           reducedTheta<=2*pi*j/numBins);
        reducedDistMed(j,i) = median(reducedVector(checkInd));
        reducedDistMad(j,i) = mad(reducedVector(checkInd));
        
        checkInd = find(reducedCorrected>=2*pi*(j-1)/numBins & ...
           reducedCorrected<=2*pi*j/numBins);
        reducedCorrDistMed(j,i) = median(reducedVector(checkInd));
        reducedCorrDistMad(j,i) = mad(reducedVector(checkInd));
    end
    thetaBins = [2*pi/numBins:2*pi/numBins:2*pi] - 2*pi/numBins/2;
    
    
    thetaMedian(i) = median(theta);
    thetaMad(i) = mad(theta);
    correctedMedian(i) = median(correctedTheta);
    correctedMad(i) = mad(correctedTheta);
    checkVec = reshape(reducedTheta,numel(reducedTheta),1);
    checkVec = checkVec(isnan(checkVec)==0);
    reducedMedian(i) = median(checkVec);
    reducedMad(i) = mad(checkVec);
    checkVec = reshape(reducedCorrected,numel(reducedCorrected),1);
    checkVec = checkVec(isnan(checkVec)==0);
    reducedCorMedian(i) = median(checkVec);
    reducedCorMad(i) = mad(checkVec);
    
    
    figure(1); clf;
    set(gcf,'Units','normalized','OuterPosition',[0.1 0.07 0.7 0.9]);
    subplot(2,2,1);
        plot(theta*180/pi,vector,'b.');
        hold on;
        errorbar(thetaBins*180/pi,vectorDistMed(:,i),vectorDistMad(:,i),...
            'r.-','Linewidth',2);
        ylabel('v (m/s)'); xlabel('\theta (deg)');
    subplot(2,2,2);
        plot(correctedTheta*180/pi,vector,'b.');
        hold on;
        errorbar(thetaBins*180/pi,correctedDistMed(:,i),correctedDistMad(:,i),...
            'r.-','Linewidth',2);
        ylabel('v (m/s)'); xlabel('\theta (deg)');
    subplot(2,2,3);
        plot(reducedTheta*180/pi,reducedVector,'b.');
        hold on;
        errorbar(thetaBins*180/pi,reducedDistMed(:,i),reducedDistMad(:,i),...
            'r.-','Linewidth',2);
        ylabel('v (m/s)'); xlabel('\theta (deg)');
    subplot(2,2,4);
        plot(reducedCorrected*180/pi,reducedVector,'b.');
        hold on;
        errorbar(thetaBins*180/pi,reducedCorrDistMed(:,i),reducedCorrDistMad(:,i),...
            'r.-','Linewidth',2);
        ylabel('v (m/s)'); xlabel('\theta (deg)');
    if saveArg == 1
        print(gcf,sprintf('OutputV2\\ThetaChart_%s',listDir(i).name),...
            '-djpeg','-r300');
    end
    
    figure(2); clf;
    set(gcf,'Units','normalized','OuterPosition',[0.15 0.1 0.65 0.85]);
    subplot(2,2,1);
        hist(theta,72); hold on;
        plot(thetaMean(i)*[1,1],get(gca,'YLim'),'k-');
        plot((thetaMean(i)-thetaDev(i))*[1,1],get(gca,'YLim'),'k--');
        plot((thetaMean(i)+thetaDev(i))*[1,1],get(gca,'YLim'),'k--');
        plot(thetaMedian(i)*[1,1],get(gca,'YLim'),'r-');
        plot((thetaMedian(i)-thetaMad(i))*[1,1],get(gca,'YLim'),'r--');
        plot((thetaMedian(i)+thetaMad(i))*[1,1],get(gca,'YLim'),'r--');
    subplot(2,2,2);
        hist(correctedTheta,72); hold on;
        plot(correctedMean(i)*[1,1],get(gca,'YLim'),'k-');
        plot((correctedMean(i)-correctedDev(i))*[1,1],get(gca,'YLim'),'k--');
        plot((correctedMean(i)+correctedDev(i))*[1,1],get(gca,'YLim'),'k--');
        plot(correctedMedian(i)*[1,1],get(gca,'YLim'),'r-');
        plot((correctedMedian(i)-correctedMad(i))*[1,1],get(gca,'YLim'),'r--');
        plot((correctedMedian(i)+correctedMad(i))*[1,1],get(gca,'YLim'),'r--');
    subplot(2,2,3);
        hist(reshape(reducedTheta,numel(reducedTheta),1),72); hold on;
        plot(reducedMean(i)*[1,1],get(gca,'YLim'),'k-');
        plot((reducedMean(i)-reducedDev(i))*[1,1],get(gca,'YLim'),'k--');
        plot((reducedMean(i)+reducedDev(i))*[1,1],get(gca,'YLim'),'k--');
        plot(reducedMedian(i)*[1,1],get(gca,'YLim'),'r-');
        plot((reducedMedian(i)-reducedMad(i))*[1,1],get(gca,'YLim'),'r--');
        plot((reducedMedian(i)+reducedMad(i))*[1,1],get(gca,'YLim'),'r--');
    subplot(2,2,4);
        hist(reshape(reducedCorrected,numel(reducedCorrected),1),72); hold on;
        plot(reducedCorMean(i)*[1,1],get(gca,'YLim'),'k-');
        plot((reducedCorMean(i)-reducedCorDev(i))*[1,1],get(gca,'YLim'),'k--');
        plot((reducedCorMean(i)+reducedCorDev(i))*[1,1],get(gca,'YLim'),'k--');
        plot(reducedCorMedian(i)*[1,1],get(gca,'YLim'),'r-');
        plot((reducedCorMedian(i)-reducedCorMad(i))*[1,1],get(gca,'YLim'),'r--');
        plot((reducedCorMedian(i)+reducedCorMad(i))*[1,1],get(gca,'YLim'),'r--');
    if saveArg == 1
        print(gcf,sprintf('OutputV2\\ThetaHist_%s',listDir(i).name),...
            '-djpeg','-r300');
    end
    
    pause(0.1);
end


%% Output

figure(3); clf;
set(gcf,'Units','normalized','OuterPosition',[0.05 0.15 0.85 0.55]);
subplot(1,3,1);
    errorbar(thetaMedian*180/pi,thetaMad*180/pi,'b.-','Markersize',10); hold on;
    errorbar(correctedMedian*180/pi,correctedMad*180/pi,'g.-','Markersize',10);
    errorbar(reducedMedian*180/pi,reducedMad*180/pi,'r.-','Markersize',10);
    errorbar(reducedCorMedian*180/pi,reducedCorMad*180/pi,'k.-','Markersize',10);
    xlabel('Frame #'); ylabel('\theta (deg)'); title('Mean, Deviation');
subplot(1,3,2);
    plot(thetaMedian*180/pi,'b.-','Markersize',10); hold on;
    plot(correctedMedian*180/pi,'g.-','Markersize',10);
    plot(reducedMedian*180/pi,'r.-','Markersize',10);
    plot(reducedCorMedian*180/pi,'k.-','Markersize',10);
    xlabel('Frame #'); title('Mean');
subplot(1,3,3);
    plot(thetaMad*180/pi,'b.-','Markersize',10); hold on;
    plot(correctedMad*180/pi,'g.-','Markersize',10);
    plot(reducedMad*180/pi,'r.-','Markersize',10);
    plot(reducedCorMad*180/pi,'k.-','Markersize',10);
    xlabel('Frame #'); title('Deviation');
if saveArg == 1
    print(gcf,'Theta_Dist','-djpeg','-r300');
end

if saveArg == 1
    save('Output.mat','time','B','L','x','y','vx','vy',...
        'corrected*','reduced*','theta*','vector*','intrudeErupt','velocityMean');
end



